I use the spaero::get_stats() function to calculate EWS over a moving window for the measles case incidence data from four cities in Niger. I use a bandwidth of 35 weeks, resulting in a 70 week window. All EWS are calculated with backward_only = TRUE so that EWS are only calculated based on data in the past.

# Load packages -----------------------------------------------------------
library(tidyverse)
library(spaero)
# Load data ---------------------------------------------------------------
fname <- "../../data/clean-data/weekly-measles-incidence-niger-cities-clean.RDS"
measles_data <- readRDS(fname) %>%
  filter(year > 1994)  # drop first NA year, only used for modeling
# Calculate EWS -----------------------------------------------------------
all_stats <- tibble()  # empty tibble for storage
for(do_region in unique(measles_data$region)){
  
  cases <- measles_data %>%
    filter(region == do_region) %>%
    pull(cases)
  
  city_stats <- spaero::get_stats(
    x = cases,
    center_trend = "local_constant", 
    center_kernel = "uniform", 
    center_bandwidth = 35, 
    stat_trend = "local_constant", 
    stat_kernel = "uniform", 
    stat_bandwidth = 35, 
    lag = 1, 
    backward_only = TRUE
  )$stats
  
  city_stats_tb <- as_tibble(city_stats) %>%
    mutate(
      time_iter = 1:n(),
      date = unique(measles_data$date)
    ) %>%
    gather(key = ews, value = value, -time_iter, -date) %>%
    mutate(region = do_region)
  
  all_stats <- bind_rows(all_stats, city_stats_tb)
}

Define a function to plot the EWS and the observations.

plot_it <- function(dates, cases, ews, lab, title){
  par(mar = c(5,5,2,5))
  plot(dates, cases, type = "l", main = title)
  par(new = T)
  plot(dates, ews, type = "l", col = "red", axes=F, xlab=NA, ylab=NA)
  axis(side = 4)
  mtext(side = 4, line = 3, lab, cex = 0.7)
}

Agadez

do_region <- "Agadez (City)"
dates <- unique(measles_data$date) 
par(mfrow = c(4,3))
for(do_ews in unique(all_stats$ews)){
  cases <- filter(measles_data, region == do_region) %>% pull(cases)
  ews <- filter(all_stats, region == do_region & ews == do_ews) %>% pull(value)
  plot_it(dates, cases, ews, lab = "ews", title = do_ews)
}

Maradi

do_region <- "Maradi (City)"
dates <- unique(measles_data$date) 
par(mfrow = c(4,3))
for(do_ews in unique(all_stats$ews)){
  cases <- filter(measles_data, region == do_region) %>% pull(cases)
  ews <- filter(all_stats, region == do_region & ews == do_ews) %>% pull(value)
  plot_it(dates, cases, ews, lab = "ews", title = do_ews)
}

Niamey

do_region <- "Niamey (City)"
dates <- unique(measles_data$date) 
par(mfrow = c(4,3))
for(do_ews in unique(all_stats$ews)){
  cases <- filter(measles_data, region == do_region) %>% pull(cases)
  ews <- filter(all_stats, region == do_region & ews == do_ews) %>% pull(value)
  plot_it(dates, cases, ews, lab = "ews", title = do_ews)
}

Zinder

do_region <- "Zinder (City)"
dates <- unique(measles_data$date) 
par(mfrow = c(4,3))
for(do_ews in unique(all_stats$ews)){
  cases <- filter(measles_data, region == do_region) %>% pull(cases)
  ews <- filter(all_stats, region == do_region & ews == do_ews) %>% pull(value)
  plot_it(dates, cases, ews, lab = "ews", title = do_ews)
}

Post-outbreak EWS

As the plots above show, most EWS get swamped out by the “memory” of outbreaks. To try and get around this, I am going to look at EWS after resetting them following outbreaks. I don’t think this is too contrived because real-world early detection systems would have to do the same thing.

How to ID outbreaks

Following Toby’s suggestion, I am going to look at the distribution of outbreak sizes (i.e., cumulative cases in each year).

pop_data <- tibble(
  region = c("Agadez (City)", "Maradi (City)", "Niamey (City)", "Zinder (City)"),
  pop = c(118224, 267249, 1027000, 322935)
)
cusum_data <- measles_data %>%
  group_by(region, year) %>%
  summarise(cases = sum(cases)) %>%
  left_join(pop_data)
Joining, by = "region"
ggplot(cusum_data, aes(x = cases/pop)) +
  geom_histogram(aes(fill = region))

cusum_cases <- sort(cusum_data$cases / cusum_data$pop)
fit1 <- MASS::fitdistr(cusum_cases, "exponential") 
x <- cusum_cases
mixexplik <- function(p,lambda1,lambda2) {
  z <- p*dexp(x,lambda1) + (1-p)*dexp(x,lambda2)
  return(-sum(log(z)))
}
mle_fit <- stats4::mle(minuslogl=mixexplik, 
                       start= list(p=0.2,
                                   lambda1 = as.numeric(fit1$estimate), 
                                   lambda2 = as.numeric(fit1$estimate)), 
                       method="Nelder-Mead")
NaNs producedNaNs producedNaNs produced
z1 <- rexp(1000, rate = mle_fit@coef["lambda1"])
z2 <- rexp(1000, rate = mle_fit@coef["lambda2"])
zall <- mle_fit@coef["p"]*z1 + (1-mle_fit@coef["p"])*z2
plot(density(zall))
lines(density(z1), col = "red")
lines(density(z2), col = "blue")

hist(x, breaks = 25, freq = FALSE, main = NULL, ylim = c(0,1300),
     xlab = "incidence (reports/person)", col = "grey")
lines(density(z1), col = "red", lwd = 2)
lines(density(z2), col = "blue", lwd = 2)
legend(0.0025, 800, legend=c("small outbreaks", "large outbreaks"),
       col=c("red", "blue"), lty = 1, cex = 0.8, lwd = 2)

LS0tCnRpdGxlOiAiRVdTIGluIG1lYXNsZXMgZGF0YSIKb3V0cHV0OiBodG1sX25vdGVib29rCi0tLQoKSSB1c2UgdGhlIGBzcGFlcm86OmdldF9zdGF0cygpYCBmdW5jdGlvbiB0byBjYWxjdWxhdGUgRVdTIG92ZXIgYSBtb3Zpbmcgd2luZG93IGZvciB0aGUgbWVhc2xlcyBjYXNlIGluY2lkZW5jZSBkYXRhIGZyb20gZm91ciBjaXRpZXMgaW4gTmlnZXIuCkkgdXNlIGEgYmFuZHdpZHRoIG9mIDM1IHdlZWtzLCByZXN1bHRpbmcgaW4gYSBgciAzNSoyYCB3ZWVrIHdpbmRvdy4KQWxsIEVXUyBhcmUgY2FsY3VsYXRlZCB3aXRoIGBiYWNrd2FyZF9vbmx5ID0gVFJVRWAgc28gdGhhdCBFV1MgYXJlIG9ubHkgY2FsY3VsYXRlZCBiYXNlZCBvbiBkYXRhIGluIHRoZSBwYXN0LgoKYGBge3IgY2FsYy1ld3N9CiMgTG9hZCBwYWNrYWdlcyAtLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLQoKbGlicmFyeSh0aWR5dmVyc2UpCmxpYnJhcnkoc3BhZXJvKQoKCiMgTG9hZCBkYXRhIC0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLQoKZm5hbWUgPC0gIi4uLy4uL2RhdGEvY2xlYW4tZGF0YS93ZWVrbHktbWVhc2xlcy1pbmNpZGVuY2UtbmlnZXItY2l0aWVzLWNsZWFuLlJEUyIKbWVhc2xlc19kYXRhIDwtIHJlYWRSRFMoZm5hbWUpICU+JQogIGZpbHRlcih5ZWFyID4gMTk5NCkgICMgZHJvcCBmaXJzdCBOQSB5ZWFyLCBvbmx5IHVzZWQgZm9yIG1vZGVsaW5nCgoKIyBDYWxjdWxhdGUgRVdTIC0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tCgphbGxfc3RhdHMgPC0gdGliYmxlKCkgICMgZW1wdHkgdGliYmxlIGZvciBzdG9yYWdlCgpmb3IoZG9fcmVnaW9uIGluIHVuaXF1ZShtZWFzbGVzX2RhdGEkcmVnaW9uKSl7CiAgCiAgY2FzZXMgPC0gbWVhc2xlc19kYXRhICU+JQogICAgZmlsdGVyKHJlZ2lvbiA9PSBkb19yZWdpb24pICU+JQogICAgcHVsbChjYXNlcykKICAKICBjaXR5X3N0YXRzIDwtIHNwYWVybzo6Z2V0X3N0YXRzKAogICAgeCA9IGNhc2VzLAogICAgY2VudGVyX3RyZW5kID0gImxvY2FsX2NvbnN0YW50IiwgCiAgICBjZW50ZXJfa2VybmVsID0gInVuaWZvcm0iLCAKICAgIGNlbnRlcl9iYW5kd2lkdGggPSAzNSwgCiAgICBzdGF0X3RyZW5kID0gImxvY2FsX2NvbnN0YW50IiwgCiAgICBzdGF0X2tlcm5lbCA9ICJ1bmlmb3JtIiwgCiAgICBzdGF0X2JhbmR3aWR0aCA9IDM1LCAKICAgIGxhZyA9IDEsIAogICAgYmFja3dhcmRfb25seSA9IFRSVUUKICApJHN0YXRzCiAgCiAgY2l0eV9zdGF0c190YiA8LSBhc190aWJibGUoY2l0eV9zdGF0cykgJT4lCiAgICBtdXRhdGUoCiAgICAgIHRpbWVfaXRlciA9IDE6bigpLAogICAgICBkYXRlID0gdW5pcXVlKG1lYXNsZXNfZGF0YSRkYXRlKQogICAgKSAlPiUKICAgIGdhdGhlcihrZXkgPSBld3MsIHZhbHVlID0gdmFsdWUsIC10aW1lX2l0ZXIsIC1kYXRlKSAlPiUKICAgIG11dGF0ZShyZWdpb24gPSBkb19yZWdpb24pCiAgCiAgYWxsX3N0YXRzIDwtIGJpbmRfcm93cyhhbGxfc3RhdHMsIGNpdHlfc3RhdHNfdGIpCn0KCmBgYAoKRGVmaW5lIGEgZnVuY3Rpb24gdG8gcGxvdCB0aGUgRVdTIGFuZCB0aGUgb2JzZXJ2YXRpb25zLgoKYGBge3IgcGxvdC1ld3MtZnVuY30KcGxvdF9pdCA8LSBmdW5jdGlvbihkYXRlcywgY2FzZXMsIGV3cywgbGFiLCB0aXRsZSl7CiAgcGFyKG1hciA9IGMoNSw1LDIsNSkpCiAgcGxvdChkYXRlcywgY2FzZXMsIHR5cGUgPSAibCIsIG1haW4gPSB0aXRsZSkKICBwYXIobmV3ID0gVCkKICBwbG90KGRhdGVzLCBld3MsIHR5cGUgPSAibCIsIGNvbCA9ICJyZWQiLCBheGVzPUYsIHhsYWI9TkEsIHlsYWI9TkEpCiAgYXhpcyhzaWRlID0gNCkKICBtdGV4dChzaWRlID0gNCwgbGluZSA9IDMsIGxhYiwgY2V4ID0gMC43KQp9CmBgYAoKIyMjIEFnYWRlegoKYGBge3IgYWdhZGV6LCBmaWcuaGVpZ2h0PTh9CmRvX3JlZ2lvbiA8LSAiQWdhZGV6IChDaXR5KSIKZGF0ZXMgPC0gdW5pcXVlKG1lYXNsZXNfZGF0YSRkYXRlKSAKcGFyKG1mcm93ID0gYyg0LDMpKQpmb3IoZG9fZXdzIGluIHVuaXF1ZShhbGxfc3RhdHMkZXdzKSl7CiAgY2FzZXMgPC0gZmlsdGVyKG1lYXNsZXNfZGF0YSwgcmVnaW9uID09IGRvX3JlZ2lvbikgJT4lIHB1bGwoY2FzZXMpCiAgZXdzIDwtIGZpbHRlcihhbGxfc3RhdHMsIHJlZ2lvbiA9PSBkb19yZWdpb24gJiBld3MgPT0gZG9fZXdzKSAlPiUgcHVsbCh2YWx1ZSkKICBwbG90X2l0KGRhdGVzLCBjYXNlcywgZXdzLCBsYWIgPSAiZXdzIiwgdGl0bGUgPSBkb19ld3MpCn0KYGBgCgojIyMgTWFyYWRpCmBgYHtyIG1hcmFkaSwgZmlnLmhlaWdodD04fQpkb19yZWdpb24gPC0gIk1hcmFkaSAoQ2l0eSkiCmRhdGVzIDwtIHVuaXF1ZShtZWFzbGVzX2RhdGEkZGF0ZSkgCnBhcihtZnJvdyA9IGMoNCwzKSkKZm9yKGRvX2V3cyBpbiB1bmlxdWUoYWxsX3N0YXRzJGV3cykpewogIGNhc2VzIDwtIGZpbHRlcihtZWFzbGVzX2RhdGEsIHJlZ2lvbiA9PSBkb19yZWdpb24pICU+JSBwdWxsKGNhc2VzKQogIGV3cyA8LSBmaWx0ZXIoYWxsX3N0YXRzLCByZWdpb24gPT0gZG9fcmVnaW9uICYgZXdzID09IGRvX2V3cykgJT4lIHB1bGwodmFsdWUpCiAgcGxvdF9pdChkYXRlcywgY2FzZXMsIGV3cywgbGFiID0gImV3cyIsIHRpdGxlID0gZG9fZXdzKQp9CmBgYAoKIyMjIE5pYW1leQpgYGB7ciBuaWFtZXksIGZpZy5oZWlnaHQ9OH0KZG9fcmVnaW9uIDwtICJOaWFtZXkgKENpdHkpIgpkYXRlcyA8LSB1bmlxdWUobWVhc2xlc19kYXRhJGRhdGUpIApwYXIobWZyb3cgPSBjKDQsMykpCmZvcihkb19ld3MgaW4gdW5pcXVlKGFsbF9zdGF0cyRld3MpKXsKICBjYXNlcyA8LSBmaWx0ZXIobWVhc2xlc19kYXRhLCByZWdpb24gPT0gZG9fcmVnaW9uKSAlPiUgcHVsbChjYXNlcykKICBld3MgPC0gZmlsdGVyKGFsbF9zdGF0cywgcmVnaW9uID09IGRvX3JlZ2lvbiAmIGV3cyA9PSBkb19ld3MpICU+JSBwdWxsKHZhbHVlKQogIHBsb3RfaXQoZGF0ZXMsIGNhc2VzLCBld3MsIGxhYiA9ICJld3MiLCB0aXRsZSA9IGRvX2V3cykKfQpgYGAKCiMjIyBaaW5kZXIKYGBge3IgemluZGVyLCBmaWcuaGVpZ2h0PTh9CmRvX3JlZ2lvbiA8LSAiWmluZGVyIChDaXR5KSIKZGF0ZXMgPC0gdW5pcXVlKG1lYXNsZXNfZGF0YSRkYXRlKSAKcGFyKG1mcm93ID0gYyg0LDMpKQpmb3IoZG9fZXdzIGluIHVuaXF1ZShhbGxfc3RhdHMkZXdzKSl7CiAgY2FzZXMgPC0gZmlsdGVyKG1lYXNsZXNfZGF0YSwgcmVnaW9uID09IGRvX3JlZ2lvbikgJT4lIHB1bGwoY2FzZXMpCiAgZXdzIDwtIGZpbHRlcihhbGxfc3RhdHMsIHJlZ2lvbiA9PSBkb19yZWdpb24gJiBld3MgPT0gZG9fZXdzKSAlPiUgcHVsbCh2YWx1ZSkKICBwbG90X2l0KGRhdGVzLCBjYXNlcywgZXdzLCBsYWIgPSAiZXdzIiwgdGl0bGUgPSBkb19ld3MpCn0KYGBgCgojIFBvc3Qtb3V0YnJlYWsgRVdTCgpBcyB0aGUgcGxvdHMgYWJvdmUgc2hvdywgbW9zdCBFV1MgZ2V0IHN3YW1wZWQgb3V0IGJ5IHRoZSAibWVtb3J5IiBvZiBvdXRicmVha3MuClRvIHRyeSBhbmQgZ2V0IGFyb3VuZCB0aGlzLCBJIGFtIGdvaW5nIHRvIGxvb2sgYXQgRVdTIGFmdGVyIHJlc2V0dGluZyB0aGVtIGZvbGxvd2luZyBvdXRicmVha3MuCkkgZG9uJ3QgdGhpbmsgdGhpcyBpcyB0b28gY29udHJpdmVkIGJlY2F1c2UgcmVhbC13b3JsZCBlYXJseSBkZXRlY3Rpb24gc3lzdGVtcyB3b3VsZCBoYXZlIHRvIGRvIHRoZSBzYW1lIHRoaW5nLgoKIyMgSG93IHRvIElEIG91dGJyZWFrcwpGb2xsb3dpbmcgVG9ieSdzIHN1Z2dlc3Rpb24sIEkgYW0gZ29pbmcgdG8gbG9vayBhdCB0aGUgZGlzdHJpYnV0aW9uIG9mIG91dGJyZWFrIHNpemVzIChpLmUuLCBjdW11bGF0aXZlIGNhc2VzIGluIGVhY2ggeWVhcikuCgpgYGB7cn0KcG9wX2RhdGEgPC0gdGliYmxlKAogIHJlZ2lvbiA9IGMoIkFnYWRleiAoQ2l0eSkiLCAiTWFyYWRpIChDaXR5KSIsICJOaWFtZXkgKENpdHkpIiwgIlppbmRlciAoQ2l0eSkiKSwKICBwb3AgPSBjKDExODIyNCwgMjY3MjQ5LCAxMDI3MDAwLCAzMjI5MzUpCikKCmN1c3VtX2RhdGEgPC0gbWVhc2xlc19kYXRhICU+JQogIGdyb3VwX2J5KHJlZ2lvbiwgeWVhcikgJT4lCiAgc3VtbWFyaXNlKGNhc2VzID0gc3VtKGNhc2VzKSkgJT4lCiAgbGVmdF9qb2luKHBvcF9kYXRhKQoKZ2dwbG90KGN1c3VtX2RhdGEsIGFlcyh4ID0gY2FzZXMvcG9wKSkgKwogIGdlb21faGlzdG9ncmFtKGFlcyhmaWxsID0gcmVnaW9uKSkKCmN1c3VtX2Nhc2VzIDwtIHNvcnQoY3VzdW1fZGF0YSRjYXNlcyAvIGN1c3VtX2RhdGEkcG9wKQpmaXQxIDwtIE1BU1M6OmZpdGRpc3RyKGN1c3VtX2Nhc2VzLCAiZXhwb25lbnRpYWwiKSAKCnggPC0gY3VzdW1fY2FzZXMKbWl4ZXhwbGlrIDwtIGZ1bmN0aW9uKHAsbGFtYmRhMSxsYW1iZGEyKSB7CiAgeiA8LSBwKmRleHAoeCxsYW1iZGExKSArICgxLXApKmRleHAoeCxsYW1iZGEyKQogIHJldHVybigtc3VtKGxvZyh6KSkpCn0KCm1sZV9maXQgPC0gc3RhdHM0OjptbGUobWludXNsb2dsPW1peGV4cGxpaywgCiAgICAgICAgICAgICAgICAgICAgICAgc3RhcnQ9IGxpc3QocD0wLjIsCiAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgbGFtYmRhMSA9IGFzLm51bWVyaWMoZml0MSRlc3RpbWF0ZSksIAogICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgIGxhbWJkYTIgPSBhcy5udW1lcmljKGZpdDEkZXN0aW1hdGUpKSwgCiAgICAgICAgICAgICAgICAgICAgICAgbWV0aG9kPSJOZWxkZXItTWVhZCIpCgp6MSA8LSByZXhwKDEwMDAsIHJhdGUgPSBtbGVfZml0QGNvZWZbImxhbWJkYTEiXSkKejIgPC0gcmV4cCgxMDAwLCByYXRlID0gbWxlX2ZpdEBjb2VmWyJsYW1iZGEyIl0pCnphbGwgPC0gbWxlX2ZpdEBjb2VmWyJwIl0qejEgKyAoMS1tbGVfZml0QGNvZWZbInAiXSkqejIKcGxvdChkZW5zaXR5KHphbGwpKQpsaW5lcyhkZW5zaXR5KHoxKSwgY29sID0gInJlZCIpCmxpbmVzKGRlbnNpdHkoejIpLCBjb2wgPSAiYmx1ZSIpCgpoaXN0KHgsIGJyZWFrcyA9IDI1LCBmcmVxID0gRkFMU0UsIG1haW4gPSBOVUxMLCB5bGltID0gYygwLDEzMDApLAogICAgIHhsYWIgPSAiaW5jaWRlbmNlIChyZXBvcnRzL3BlcnNvbikiLCBjb2wgPSAiZ3JleSIpCmxpbmVzKGRlbnNpdHkoejEpLCBjb2wgPSAicmVkIiwgbHdkID0gMikKbGluZXMoZGVuc2l0eSh6MiksIGNvbCA9ICJibHVlIiwgbHdkID0gMikKbGVnZW5kKDAuMDAyNSwgODAwLCBsZWdlbmQ9Yygic21hbGwgb3V0YnJlYWtzIiwgImxhcmdlIG91dGJyZWFrcyIpLAogICAgICAgY29sPWMoInJlZCIsICJibHVlIiksIGx0eSA9IDEsIGNleCA9IDAuOCwgbHdkID0gMikKCmBgYAoKCg==